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EXECUTIVE  SUMMARY 

During  the  Phase  I  research  effort,  considerable  progress  was  made  in  atmospheric 
temperature  profile  determination  by  extension  of  differential  inversion  (Dl)  techniques  and  the 
application  of  the  Optical  Measure  Theory  (OMT).  These  advances  include  the  first-time-ever 
determination  of  partial  atmospheric  temperature  profiles  by  application  of  the  standard  DI 
theory  (first  applied  successfully  to  a  single  altitude  by  workers  at  Creative  Optics,  Inc.),  and 
the  first  bias-free  determinations  of  atmospheric  temperature  profiles  by  the  application  of 
OMT. 


Phase  I  research  demonstrated  the  applicability  of  the  Nonlinear  Hyperbolic  Algorithm 
(NHA)  and  OMT  to  atmospheric  profile  determination  and  produced  a  successful  temperature 
profile  using  these  methods  applied  to  TOVS  data.  See  Figure  1.  NHA/OMT  are  now  well- 
characterized  with  regard  to  the  numerical  analytic  properties  relevant  to  practical 
implementation.  As  a  result  we  now  have  a  good  theoretical  understanding  of  the  classic 
instability  problems  inherent  to  the  inversion  of  upwelling  radiances  to  obtain  atmospheric 
temperature  profiles.  These  instabilities  are  now  localized  to  a  few  algorithmic  steps  within  the 
NHA  and  are  amenable  to  control  by  well  understood  numerical  techniques.  The  behavior  of 
the  NHA  in  the  presence  of  noise  has  been  studied  with  a  view  ultimately  to  characterizing  the 
information  content  of  radiance  data  and  the  corresponding  atmospheric  structure  parameters. 

Determinations  of  partial  atmospheric  temperatures  (three  altitude  points  in  the 
neighborhood  of  the  tropopause)  were  also  obtained  using  standard  DI  techniques  by  truncating 
the  DI  series  at  appropriate  order.  These  partial  profiles  are  of  great  interest  in  that  they  seem 
to  replicate,  faithfully,  the  actual  temperature  structure  near  the  tropopause;  this  information  is 
of  great  operational  significance.  See  Figure  2. 

The  result  of  greatest  immediate  practical  importance  is  the  first  successful 
determination  of  atmospheric  temperature  profiles  using  NHA/OMT.  While  still  preliminary, 
these  results  are  of  great  significance  because: 

•  They  are  the  first  bias-free  determinations  of  the  earth’s  atmospheric  temperature 

profile. 

•  They  exhibit  a  meteorologic  character;  in  particular,  they  correctly  place  the  height 

of  the  tropopause  when  compared  with  radiosonde  data. 

•  They  demonstrate  the  theoretical  validity  of  the  NHA/OMT  approach  to 

atmospheric  temperature  profiling. 

In  addition  to  our  successful  utilization  of  TOVS  data,  two  other  two  potential  data  sets 
for  future  analysis  by  DI  and  NHA/OMT  techniques  were  analyzed  and  evaluated.  Data  from 
the  Defense  Meteorological  Satellite  Program  (DMSP)  microwave  radiometer  is  reduceable  using 
the  same  techniques  we  pioneered  for  TOVS,  once  data  conversion  and  formatting  problems  are 
solved.  An  approach  also  has  been  developed  for  initial  analysis  of  data  from  the  High- 
Resolution  Interferometer  Spectrometer  (HIS).  HIS  data  poses  new  questions,  and  an 
opportunity  to  develop  DI  and  NHA/OMT  for  new  applications. 
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Figure  1.  Bias-Free  Inverse  Tansformation.  Researchers  at  Creative  Optics,  Inc.  have  produced 
the  first  bias-free  atmospheric  temperature  profiles  from  satellite  radiance  data  by  an 
application  of  Dr.  J,  I.  F.  King's  Optical  Measure  Theory  (OMT).  The  profile  shown  here  is 
obtained  from  data  from  the  TIROS  Operational  Vertical  Sounder  (TOVS)  which  yielded  six 
atmospheric  structure  parameters.  The  general  meteorological  character  of  the  temperature 
profile  is  evident.  High  spectral  resolution  sounder  data  will  permit  calculation  of  more 
structure  parameters  to  represent  detailed  features  of  the  temperature  profile. 
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1.1  Background.  Satellite  observations  provide  the  only  practical  means  of  obtaining 
temperature  profiles  in  the  atmosphere  over  large  geographic  areas  (especially  over  oceans)  and 
over  long  periods  of  time.  This  implies  that  the  development  of  new  and  computationally 
effective  algorithms  for  measuring  atmospheric  temperature  profiles  from  satellite  data  is  a 
significant  and  important  task. 

The  measurable  quantity  most  accessible  to  satellite  observations  is  the  infrared 
upwelling  radiance.  Most  typically  this  is  measured  at  wavelength  -  15  pm  near  the 
atmospheric  C02  absorption  band.  Different  spectral  regions  in  this  wavelength  region  become 
optically  thick  at  different  heights  in  the  atmosphere,  and  so  the  upwelling  infrared  radiance  in 
those  wavelength  intervals  contains,  at  least  in  principle,  information  about  the  Planck  function 
at  different  heights  in  the  atmosphere.  On  this  basis,  one  should  be  able  to  determine  a 
temperature  profile  of  the  atmosphere  using  measurements  of  the  upwelling  infrared  radiance. 

In  practice,  we  must  solve  the  radiative  transfer  equation  for  the  Planck  function  over 
this  spectral  interval,  a  Fredholm  integral  equation  of  the  first  kind.  Since  we  are  dealing  with 
an  integral  equation,  we  must  cope  with  problems  of  non-uniqueness  of  the  solution  obtained. 
This  is  a  problem  of  practical  importance  because  the  range  of  plausible  solutions  of  the 
transfer  equation  (within  limits  imposed  by  the  accuracy  of  the  experimental  data)  cover  an 
unacceptably  large  range  of  temperatures. 

A  number  of  workers  have  approached  this  problem  with  techniques  designed  to  obtain 
the  “true”  temperature  profile  (or  at  least  a  pragmatically  useful  profile)  by  a  variety  of 
methods,  most  of  which  require  some  a  priori  assumptions  such  as  an  initial  temperature  profile. 
Application  of  relaxation  techniques  (Chahine  1970,  1972;  Smith  1970;Yeh,  Vonder  Haar,  and 
Liou  1985)  requires  the  selection  of  an  initial  temperature  profile  to  commence  an  iterative 
process.  Application  of  the  method  of  Backus  and  Gilbert  (Conrath  1972)  requires  the 
assumption  of  a  reference  temperature  profile  in  order  to  evaluate  radiative  transfer  kernels. 

Other  relevant  work  has  been  undertaken  (Twomey  1970)  to  delineate  the  limits  in 
information  content  of  solutions  of  the  radiative  transfer  equation.  In  this  vein,  Fourier  analysis 
of  the  radiative  transfer  equation  (Gautier  and  Revah  1975)  allows  treatment  of  the  temperature 
profile  as  a  band- limited  spectrum. 

In  this  report  we  consider  new  approaches  to  this  problem  based  on  Laplace  transform 
theory  (King,  Hohlfeld  and  Killian  (1988);  King  1983,  1985;  Leon  and  King  1988)— Differential 
Inversion  (DI)  and  Optical  Measure  Theory  (OMT).  The  techniques  use  the  upwelling  radiance 
measurements  and  their  logarithmic  pressure  derivatives  to  construct  an  atmospheric  temperature 
profile.  The  methods  have  the  particular  theoretical  attraction  of  assuming  no  a  priori 
information  about  the  temperature  profile  beyond  the  fact  that  it  is  reasonably  smooth  ( i.e . 
analytic).  Our  numerical  studies  and  simulations  using  the  algorithms  have  investigated  the 
viability  of  these  approaches  and  indicated  some  of  their  sensitivities  to  truncation  error, 
roundoff  error,  experimental  uncertainty,  and  other  noise  processes.  DI  and  OMT  have  been 
successfully  applied  to  a  sample  of  data  from  the  NOAA  TIROS  Operational  Vertical  Sounder 
(TOVS). 

1.2  Phase  1  Technical  Objectives.  Broadly  speaking  the  work  undertaken  during  the  Phase  I 
research  effort  was  directed  toward  adapting  existing  software  for  determination  of  atmospheric 
temperature  profiles  by  DI  to  the  specific  application  of  the  reduction  of  radiance  data.  Not 
only  did  we  achieve  our  original  objectives,  we  achieved  successes  in  new  areas  not  originally 
proposed,  viz.,  new  software  developed  for  the  implementation  of  the  OMT  and  the  Nonlinear 


Hyperbolic  Algorithm  (NHA)  which  permits  computation  of  bias-free  atmospheric  temperature 
profiles  from  radiance  data. 
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INVERSION 


The  differential  inversion  algorithm  of  King  (1988)  has  been  applied  by  research 
workers  at  Creative  Optics,  Inc.  to  obtain  the  first  bias-free  determinations  of  atmospheric 
temperature  (King,  Hohlfeld,  s»d  Kilian  1988).  The  first  such  application  used  seven  channels 
of  4.3  fim  infrared  data  from  the  Tiros  Operational  Vertical  Sounder  (TOVS)  instrument  to 
obtain  a  single  temperature  value  at  398  millibars. 

Temperature  values  in  the  differential  inversion  theory  are  expressed  in  terms  of  a  series 
expansion  involving  the  radiance  and  radiance  derivatives, 

B(z)  -  f  Ak  D(k)  R(D  (I) 

k=0 

where  B  is  the  Planck  function,  R  is  the  radiance,  z  -  -  In  P,  where  P  is  the  pressure  at  which 
the  weight  function  reaches  its  maximum  value,  and  ?  ■  -  In  p  is  a  height- like  variable  where  p 
is  barometric  pressure.  The  A('s  are  the  differential  inversion  coefficients  determined  as 
moments  of  the  weight  functions.  Originally,  all  seven  TOVS  channels  were  employed  to  obtain 
radiance  numerical  derivatives  via  centered  divided  differences  out  to  the  third  derivative  (i.e. 
the  term  containing  A4),  which  permitted  the  evaluation  of  a  single  atmospheric  temperature 
value  at  the  central  P  value  among  the  seven  channels,  as  shown  in  Figure  2. 

It  was  noted,  however,  that  convergence  of  this  series  was  very  rapid,  and  the 
contribution  of  the  highest  order  terms  was  less  than  the  expected  errors  arising  due  to 
numerical  and  experimental  sources  of  error.  Therefore,  we  truncated  the  differential  inversion 
series  at  second  order,  so  that  two  additional  atmospheric  temperature  values  could  be  obtained, 
while  still  computing  all  numerical  derivatives  as  centered  divided  differences  (Hohlfeld,  Kilian, 
Drueding,  and  Ebersole  1988).  These  additional  temperature  points  were  obtained  at  pressure 
values  of  174  and  661  millibars,  as  shown  in  Figure  3. 

A  collection  of  Zone  1  (latitude  -  60*  N  to  90*  N)  TOVS  data  were  analyzed  in  this 
fashion  and  the  results  are  collected  in  Appendix  A.  It  may  be  seen  by  comparison  with  the 
RAOB  data,  that  significant  atmospheric  structure  has  been  represented  by  these  preliminary 
results.  Most  notably,  data  points  on  either  side  of  the  tropopause  have  been  obtained, 
demonstrating  the  potential  utility  of  this  analysis  technique  for  satellite  remote  measurement  of 
the  altitude  of  the  tropopause. 
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DIFFERENTIAL  INVERSION  APPLIED  TO  TOVS  DATA  TO  OBTAIN  A 

THREE  POINT  PROFILE 


Pressure(mb) 


Temperature(K) 


Figure  3.  Differential  Inversion  Applied  to  TOVS  Data  to  Obtain  a  Three-Point  Temperature 
Profile.  Shown  here  is  a  representative  example  of  Creative  Optics’  results  for  partial 
atmospheric  temperature  profiles  (at  three  altitudes)  from  4.3  /rm  to  TOVS  data.  Note  that 
temperature  values  on  either  side  of  the  tropopause  have  been  obtained.  This  demonstrates  the 
potential  ability  of  this  algorithm  for  the  important  operational  problem  of  determination  of 
tropopause  altitude  by  remote  sensing.  Other  examples  are  presented  in  Appendix  A. 


3.1  Introduction.  Optical  Measure  Theory  (OMT)  is  a  theoretical  approach  to  the  problem 
of  determination  of  atmospheric  temperatures  in  which  an  atmospheric  temperature  profile  is 
the  fundamental  theoretical  object.  That  is  to  say,  OMT  is  formulated  in  terms  of  atmospheric 
temperature  profiles  in  the  same  sense  as  DI  is  formulated  in  terms  of  obtaining  discrete 
temperature  values  in  the  atmosphere. 

OMT  has  been  developed  on  the  basis  of  advances  in  transform  theory  (Leon  and  King 
1988,  King  1988).  Specifically  it  was  noted  by  King  that  the  radiative  transfer  equation 
represents  a  generalization  of  the  Laplace  transform,  in  the  manner  that  it  relates  the  Planck 
function  profile  to  the  radiance  profile.  Based  on  this  insight,  it  is  then  natural  to  represent  the 
radiance  profile  in  some  appropriate  functional  form  (the  selection  of  the  appropriate  functional 
form  is  a  problem  discussed  below)  and  then  to  compute  the  Planck  function  profile,  and  hence 
the  atmospheric  temperature  profile,  as  an  inverse  transform  of  the  radiance  profile. 

Several  important  advantages  derive  from  this  approach.  First  of  all,  temperature 
profiles  arise  in  a  natural  manner  in  OMT.  Secondly,  if  an  intelligent  choice  is  made  regarding 
the  functional  representation  of  the  radiance  data,  these  profiles  will  necessarily  exhibit  a 
meteorological  character.  Lastly,  and  perhaps  of  the  greatest  importance,  the  classical  numerical 
instability  of  the  atmospheric  radiance  inversion  problem  is  localized  algorithmically  to  a  few 
well-understood  steps  which  may  then  be  treated  with  specialized  techniques. 

3.2  Theoretical  Development  of  Optical  Measure  Theory.  For  sake  of  completeness,  we  shall 
derive  here  the  fundamental  equations  of  OMT.  Following  King  (1988)  and  Leon  and  King 
(1988),  we  note  that  the  equation  for  upwelling  radiance 

R(M)  -  J^°B(p)W(p//i)  %  (2) 

has  the  character  of  an  integral  transform  from  B(p),  the  Planck  function  as  a  function  of 
pressure,  p,  to  the  radiance,  R(/i),  where  /j  is  a  pressure-like  transform  variable.  Here  W(p//r) 
is  the  usual  atmospheric  weight  function.  In  the  context  of  the  standard  DI  development,  n, 
would  be  viewed  as  that  pressure  value  at  which  the  weight  function  achieved  its  maximum 
value.  If  the  actual  atmospheric  weight  function  is  well-represented  as  a  generalized 
exponential  function  (King  1985;  King,  Hohlfeld,  and  Kilian  1988) 

W(p/P)  =  Wm(p/M)  -  r~ ST)  (p//i)  exp(-m(p/n)1/m],  (3) 

then  for  m  =  I,  Eq.  (2)  immediately  assumes  the  character  of  a  Laplace  transform, 

R(^)  “  ^CB(p)e'(P/'l)  dP  (4) 

i.e.  pR(p)  is  the  Laplace  transform  of  B(p)  with  respect  to  the  transform  variable  l/p.  When 
m  *  1,  the  radiance  profile  and  Planck  function  profile  are  then  related  by  a  generalization  of 
the  Laplace  transform  defined  in  terms  of  the  generalized  exponential  function.  (For  example 
with  TOVS  data  see  Figure  4.)  All  of  the  useful  analytic  properties  of  the  Laplace  transform 
are  retained  in  this  generalization  of  the  Laplace  transform  utilizing  the  generalized  exponential 
function  kernel. 


If  a  choice  of  functional  form  is  made  for  B(p)  (or  alternatively  R(/i)),  Equation  (2)  then 
immediately  implies  the  functional  form  R(^)  (alternatively  B(p)).  This  observation  is  especially 
pertinent  if  we  note  that  the  Planck  function  profile  of  a  radiative  atmosphere  is  exponential  in 
form,  i.e. 


B(p)  -  Le‘kp  (5) 

with  k  and  L  constants  (Chandrasekhar  1960).  The  corresponding  (m  =  1)  radiance  profile 
function  is  then, 


R<*0 


L 

1  +  k/i‘ 


(6) 


This  discussion  motivates  a  choice  of  the  functional  form  by  which  radiance  data  is 
represented  (King  1988;  Leon  and  King  1988), 

R(fi)  =  a  +  b/i  +  £  n^T 

where  j  hyperbolic  terms  are  included.  The  addition  of  a  linear  term  of  form  a  +  b n  (with  a 
and  b  constants)  was  found  by  Leon  and  King  to  be  of  practical  utility  in  the  representation  of 
radiance  data.  Generalization  to  the  inclusion  of  a  polynomial  of  arbitrary  order  in  n  is 
straightforward.  The  corresponding  Planck  function  profile  obtained  by  generalized  inverse 
Laplace  transformation  (for  arbitrary  m)  of  Eq.  (7)  is 

B(p)  -  a  ♦  bp  +  £  LjEJ-kjP)  W 

i=l 


where  Em(x)  is  the  generalized  exponential  function  with  parameter  m. 

The  choice  of  an  exponential  (or  generalized  exponential  form  of  the  Planck  function 
profile,  motivated  by  the  expectation  of  the  real  atmospheric  acting  at  least  in  part  as  a 
radiative  atmosphere,  indicates  that  Eq.  (8)  is  a  natural  choice  of  functional  form  for  the 
representation  of  radiance  data. 

3.3  The  Nonlinear  Hyperbolic  Algorithm.  The  theoretical  discussion  given  above  shows  the 
importance  of  fitting  upwelting  radiance  data  to  a  functional  form  as  given  in  Eq.  (7).  Leon 
and  King  (1988)  have  developed  an  algorithm,  the  Nonlinear  Hyperbolic  Algorithm  (NHA),  for 
fitting  radiance  data  to  a  formula  of  the  form  of  Eq.  (7).  As  such,  NHA  is  the  natural  starting 
point  for  any  discussion  of  the  properties  of  the  OMT.  The  NHA  allows  determination  of  the 
constants  a,  b,  and  the  set  of  L.’s  and  kj’s  in  Eq.  (7),  which  then  immediately  determines  the 
corresponding  terms  of  the  PlancV  function  profile,  given  by  Eq.  (8),  and  hence  the  atmospheric 
temperature  profile. 

We  now  outline  the  steps  of  the  NHA.  Assume  that  we  are  given  2n  =  2j  +  2  upwelling 
radiance  measurements,  which  are  to  be  fitted  to  a  formula  of  the  form  of  Eq.  (7),  i.e.,  we  wish 
to  determine  a,  b,  L.,  ....  L.,  k  ,  ....  k.,  a  total  of  2n  constants  from  2n  radiance  measurements. 
Following  Chandrasekhar  (I960),  the  radiance  is  represented  as  a  rational  function  in  n. 


R(A«) 


Jn  +  2 


di  ♦  d2M  ♦  ■  ■  ■  +  d^ji" 

*  dn+3*i  +  .  .  .  ♦  d2nnn-3  ♦  M"-1 


POO/QOO. 


(9) 
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This  should  be  viewed  as  a  linear  system  of  equations  for  the  d.’s  of  rank  2n.  Writing 
this  system  as  P(fij)  *  0(^)1^  where  /i.  is  the  pressure  of  the  weight  function  maximum  of  the 
ith  channel,  we  obtain  2n  equations  for  the  2n  unknowns,  where  the  ith  equation  is 


dt  +  d^  +  .  .  .  +dB+1M"  -  Mn+J  -  -  •  •  •-  Rjdj „^n-J  *  (10) 


This  linear  system  may  be  solved  by  any  of  the  standard  techniques  for  linear  systems  of 
equations,  such  as  singular  value  decomposition  (Press,  Flannery,  Teukolsky,  and  Vetterling 
1986).  Matching  orders  for  n°  and  /s1  in  Eqs.  (9)  and  (7),  we  obtain 

b  -  d„*i 

a  -  <*n  "  dn+1d2n  (11) 

The  remaining  coefficients  in  Eq.  (7)  may  be  determined  by  noting  that  the  locations  of 

poles  in  R(/i)  determine  values  of  the  k/s. 

kj  -  -1/Cj  (12) 

where  c^  is  the  jth  root  of  Q(j a).  The  roots  of  QOr)  may  be  found  by  any  of  the  standard  root¬ 
finding  applicable  to  polynomials,  such  as  the  Laguerre  method  (Press,  Flannery,  Teukolsky,  and 
Vetterling  1986),  which  allows  for  the  presence  of  complex  roots.  Lastly,  the  L.  values  may  be 
determined  by  solution  of  the  linear  system  of  equations 

rj  -a  (13) 

where  j  -  1,  .  .  .,  n  is  any  subset  of  the  radiance  measurements. 

3.4  Stability  of  the  Nonlinear  Hyperbolic  Algorithm.  While  the  NHA  is  well-grounded  on 
physical  principles  applicable  to  radiative  transport  in  the  earth’s  atmosphere,  practical 
difficulties  arise  in  the  numerical  stability  of  the  NHA.  Determination  of  the  coefficients  in 
Eq.  (7)  in  application  of  the  NHA  is  analogous  to  a  classic  problem  in  curve-fitting,  that  of 
fitting  exponential  curves  to  radioactive  decay  data  of  several  radionuclides,  well-known  to  be 
numerically  unstable  (Acton  1970).  This  may  be  seen  qualitatively  by  noting  that  e*x  and 
l/(l+x)  are  roughly  equal  over  the  range  0  <  x  <  1. 

These  numerical  instabilities  manifest  themselves  as  an  ill-determinedness  in  applications 
to  radiance  data  of  the  NHA  coefficients  a,  b,  kJ(  .  .  k.,  Lr  .  .  .  L..  Small  perturbations  of  the 
input  radiance  values  lead  to  large  changes  in  the  coefficients  o^  the  fit  determined  by  the 
NHA. 

It  is  instructive  to  note  that  the  transformation  of  the  radiance  series,  Eq.  (7),  to  the 
Planck  function  series,  Eq.  (8),  is  absolutely  stable.  Furthermore,  obtaining  the  temperature 
profile  from  the  Planck  function  profile  involves  a  stable  and  continuous  transformation. 
Therefore,  all  of  the  classical  numerical  instability  of  the  problem  of  determination  of  the 
atmospheric  temperature  profile  from  upwelling  atmospheric  radiances  has  been  localized 
algorithmically  in  the  NHA.  In  particular,  the  principal  unstable  operations  in  the  NHA  involve 
the  solution  of  the  linear  system  for  the  coefficients  d.,  .  .  .,  dJn.  It  is  of  great  practical 
significance  in  applications  to  inversion  problems  to  fully  understand  the  origin  of  numerical 
instabilities  in  order  to  control  them  effectively  by  specialized  numerical  techniques. 
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Stability  Properties.  Creative  Optics,  Inc.  has  carried  out  a  successful  implementation  of  the 
NHA  in  the  FORTRAN  language  running  on  PC-AT  class  microcomputers.  Previous 
implementations  by  Leon  and  King  (1988)  utilized  a  PC-based  MATHCAD  package  suitable  for 
exploratory  work,  while  the  FORTRAN-based  implementation  discussed  here  represents  a 
significant  step  towards  an  efficient,  operational  of  the  NHA. 

Several  specialized  numerical  techniques  were  employed  to  optimize  the  numerical 
stability  properties  of  the  NHA.  Most  important  is  the  use  of  singular  value  decomposition 
(Press,  Flannery,  Teukolsky,  and  Vetterling  1986)  for  the  solution  of  the  linear  system 
determining  the  coefficients  dr  .  .  .,  d2n.  Singular  value  decomposition  is  useful  for  the  solution 
of  ill-conditioned  linear  systems  because  it  is  capable  of  returning  partial  solutions  if  the  quality 
of  the  data  set  will  not  permit  a  well-conditioned  solution  of  full  rank. 

Solution  for  the  roots  of  the  polynomial  equations  for  the  k/s  was  carried  out  by 
Laguerre’s  method  (Press,  Flannery,  Teukolsky,  and  Vetterling  1986).  Laguerre’s  method  is  able 
to  compute  complex  polynomial  roots  as  effectively  as  real  roots,  and  its  convergence  to  roots  is 
guaranteed.  When  combined  with  standard  decimation-in-order  techniques,  Laguerre’s  method 
is  an  effective  approach  to  root-finding  of  polynomials,  even  polynomials  of  moderately  large 
order. 


Studies  were  undertaken  using  this  FORTRAN  implementation  of  the  NHA  to 
characterize  the  numerical  stability  of  the  NHA  in  controlled  computational  experiments.  These 
studies  involved  applying  the  NHA  both  to  synthetic  radiance  values  generated  from  integration 
of  the  radiative  transfer  equation,  supplied  with  a  realistic  atmospheric  temperature  profile,  and 
with  radiance  values  generated  by  application  of  Eq.  (7)  with  known  coefficients.  These  values 
may  be  computed  with  known  values  of  Gaussian  noise  so  as  to  determine  the  sensitivity  of  the 
NHA  to  experimental  errors. 

Studies  have  been  undertaken  to  investigate  the  sensitivity  of  the  NHA  both  to  noise  in 
the  radiance  amplitudes  and  noise  in  the  P  values,  i.e,  uncertainty  in  the  pressure  value  at  which 
the  corresponding  channel  weight  functions  achieve  their  maximum  values.  These  latter  sources 
of  error  were  motivated  by  noting  practical  difficulties  experienced  in  fitting  generalized 
exponential  functional  forms  to  some  atmospheric  weight  functions.  This  implies,  in  effect,  that 
some  P  values  are  not  well-determined  at  the  5%  level. 

During  the  Phase  I  effort  we  have  obtained  a  preliminary  quantification  of  the  level  of 
sensitivity  of  the  NHA  to  experimental  noise.  The  NHA  still  functions  reliably  and  effectively 
in  the  presence  of  Gaussian  noise  in  radiance  values  at  the  5%  to  10%  level.  It  is  interesting 
and  significant  to  note  that  the  NHA  is  comparatively  sensitive,  however,  to  noise  in  the  P 
values.  Noise  in  excess  of  1%  in  P  may  prevent  the  NHA  from  finding  physically  meaningful 
coefficients  for  the  nonlinear  hyperbolic  fit.  This  result  highlights  the  importance  of  the 
techniques  we  have  developed  at  Creative  Optics  for  the  analysis  of  atmospheric  weight 
functions,  which  can  now  be  seen  to  be  of  significance  for  the  practical  application  of  the 
NHA. 


We  have  characterized  the  behavior  of  the  NHA  under  these  circumstances  which  are 
numerically  ill-conditioned  due  to  the  presence  of  experimental  errors.  Typically  negative 
values  of  the  kj’s  are  obtained,  which  correspond  to  poles  in  the  radiance  formula  at  positive 
values  of  pressure,  or  values  of  k,’s  which  are  complex  conjugate  pairs.  Clearly,  both  of  these 
cases  correspond  to  nonphysical  solutions  and  are  readily  identified  as  such.  This  error¬ 
checking  feature  is  important  for  an  operational  algorithm  which  must  perform  reliably  in  the 
presence  of  occasional  bad  data. 


II 


Leon  and  King  (1988)  have  demonstrated  the  utility  of  the  NHA  for  detecting  spurious 
values  in  radiance  data  corrupted  by  a  single  bad  channel.  Our  studies  confirm  the  ability  of 
the  NHA  to  detect  the  corrupted  channel  and  to  generate  the  best  possible  fit  to  the  remaining 
radiance  data  when  the  bad  channel  value  is  ignored.  We  have  furthered  these  studies  uy 
considering  experimental  errors  in  which  each  channel  is  degraded  by  some  level  of  random 
noise.  Under  these  circumstances,  when  the  noise  level  in  the  radiances  reaches  the  5%  to  10% 
level  (as  discussed  above)  the  values  of  the  k/s  become  poorly  determined  (due  to  loss  of 
numerical  condition  of  the  NHA)  and  poles  are  introduced  in  apparently  random  locations  in 
pressure  space.  This  behavior  should  be  contrasted  with  the  behavior  elucidated  by  Leon  and 
King  where  the  pole  occurs  at  a  definite  location  in  proximity  to  the  single  spurious  channel. 

It  is  important  to  note  that  the  fundamental  results  of  this  study  relative  to  the 
sensitivity  of  the  NHA  to  uncertainty  in  P  values  and  the  behavior  of  the  NHA  in  the  presence 
of  noise  now  permit  practical  application  of  the  NHA  under  realistic  conditions.  This  is 
especially  true  when  the  atmospheric  weight  functions  are  well-determined  and  functionally 
well-characterized. 


12 


4.1  Motivation  and  Theoretical  Development.  Algorithmic  extension  of  the  NHA  was 
considered  as  a  means  of  coping  with  the  inherent  numerical  instability  of  the  NHA.  It  is 
desirable  to  retain  the  functional  form  of  Eq.  (7)  as  an  appropriate  representation  of  radiances 
from  a  radiative  atmosphere,  however,  this  choice  of  hyperbolic  functions  is  the  ultimate  origin 
of  the  numerical  instability  of  the  NHA.  The  numerical  instability  of  the  NHA  may  be  viewed 
as  arising  when  the  information  content  of  the  radiance  data  set  is  insufficient  to  determine  all 
the  coefficients  in  the  nonlinear  hyperbolic  representation.  A  strategy  suggested  by  this 
observation  is  to  attempt  to  fit  a  limited  set  of  hyperbolic  terms,  in  a  least-squares  sense,  to  the 
available  radiance  data. 

In  the  Phase  I  research  effort  carried  out  at  Creative  Optics,  Inc.,  we  have  implemented 
an  extension  of  the  NHA  based  on  nonlinear  least-squares  (NLLS)  techniques.  As  presently 
implemented,  the  NLLS  NHA  is  a  test-bed  algorithm  used  to  investigate  detailed  approaches  for 
developing  fits  of  the  form  of  Eq.  (7)  to  radiance  data. 

The  NLLS  NHA  is  carried  out  by  the  standard  Levenburg-Marquardt  algorithm  (Press, 
Flannery,  Teukolsky,  and  Vetterling  1986).  In  typical  operation  of  the  NLLS  NHA,  a  fit  is 
made  to  the  linear  terms  (determination  of  a  and  b)  and  a  single  hyperbolic  term  (determination 
of  Lj  and  kj).  Successive  hyperbolic  terms  are  added  until  a  nonphysical  term  is  obtained  (i.e. 
with  a  negative  k  value),  which  is  then  discarded  and  the  NLLS  NHA  algorithm  is  then 
terminated. 

4.2  Initial  Results  from  the  NLLS  NHA.  During  the  Phase  I  research  effort  at  Creative 
Optics,  Inc.,  the  NLLS  NHA  was  applied  to  a  few  TOVS  Zone  1  measurements  and  a  successful 
inversion  of  TOVS  radiance  data  to  yield  an  atmospheric  temperature  profile  was  obtained.  Six 
TOVS  channels  (1,  2,  4,  5,  6,  7)  were  used  to  obtain  the  results  shown  in  Figure  1.  This  profile 
has  an  evident  meteorological  character,  particularly  in  that  a  minimum  in  the  atmospheric 
temperature  profile  approximately  at  the  tropopause  altitude  is  seen.  It  is  also  evident  that  the 
NLLS  NHA  generates  the  smoothest  radiance  profile  consistent  with  the  observed  radiance  data. 
This  is  a  highly  desirable  feature  in  that  the  simplest  possible  atmospheric  temperature  profile  is 
constructed.  Atmospheric  structure  not  supported  by  the  radiance  data  is  not  "invented*  by  the 
algorithm.  The  six  atmospheric  structure  parameters  determined  (a,  b,  k1,  Lr  k2,  L2)  do  not 
contain  sufficient  information  to  represent  detailed  structures  in  the  atmospheric  temperature 
profile  such  as  the  ground-level  temperature  inversion  seen  in  the  radiosonde  data 
corresponding  to  this  temperature  profile  in  Figure  1.  It  is  expected  that  application  of  the 
NLLS  NHA  to  temperate  latitude  (Zone  2)  TOVS  data  will  be  more  successful  in  terms  of 
detailed  matches  between  radiosonde  data  and  temperature  profiles  derived  by  the  NLLS  NHA. 
The  preliminary  result  shown  here  demonstrates  the  potential  practical  utility  of  the  NLLS 
NHA,  particularly  for  the  important  problem  of  remote-sensing  of  the  tropopause  altitude. 

4.3  Recommended  Upgrades  to  the  NLLS  NHA.  Choice  of  Basis  Functions,  and 
Determination  of  Atmospheric  Structure  Parameters.  The  performance  of  the  present  version 
of  the  NLLS  NHA  is  limited  by  the  geometry  of  the  x 3  (error)  surface.  NLLS  algorithms 
attempt  to  determine  a  set  of  parameters  of  a  functional  fit  corresponding  to  a  (global) 
minimum  in  x3-  However,  if  the  error  surface  has  many  secondary  (local)  minima,  it  is  possible 
for  the  NLLS  algorithm  to  return  a  solution  corresponding  to  one  of  these  secondary  minima  in 
X3  rather  than  the  true  global  minimum  corresponding  to  the  desired  physical  solution.  Such 
irregular  error  surfaces  are  characteristic  of  ill-conditioned  curve  fitting  problems  (Acton  1970). 


Standard  techniques  have  been  developed  for  nonlinear  least-squares  curve  fitting  under 
such  circumstances.  These  techniques  include  grid  search,  Monte  Carlo  methods,  addition  of 
numerical  "momentum"  terms  to  NLLS,  conjugate  gradient  search  on  the  x3  surface,  and  others. 


Choice  of  an  appropriate  method  is  best  made  after  a  more  detailed  characterization  of  the 
geometric  properties  of  the  x2  surface,  which  is  appropriate  to  a  Phase  II  effort.  Application  of 
such  techniques  can  be  expected  to  make  determination  of  coefficients  of  an  NHA  fit  routine 
and  completely  independent  of  all  starting  values  used  to  initiate  the  NLLS  NHA  fit. 

Determination  of  the  NHA  parameters  can  be  viewed  as  the  determination  of  a  set  of 
atmospheric  structure  parameters  that  may  be  used  to  represent  corresponding  structure  in  the 
atmospheric  temperature  profile.  Since,  at  most,  only  as  many  structure  parameters  may  be 
obtained  as  the  number  of  channels  in  the  original  radiance  data,  it  is  clear  that  fundamental 
theoretical  limits  exist  as  to  the  amount  of  atmospheric  structure  which  can  be  represented  on 
the  basis  of  a  given  radiance  data  set.  Furthermore,  if  the  channel  measurements  of  radiance 
are  not  all  statistically  independent,  information  theoretical  principles  suggest  that  only  a  smaller 
number  of  NHA  parameters  may  be  consistently  obtained. 

These  considerations  on  the  information  content  of  radiance  data  and  their  translation 
into  atmospheric  structure  parameters  by  NHA/OMT  suggest  exciting  possibilities  for  the  design 
of  temperature  sounding  instruments  which  can  be  optimized  to  represent  structure  in  the 
atmospheric  temperature  profile  with  specified  resolution  in  height.  In  addition  to  the 
production  of  atmospheric  temperature  profiles  as  such,  this  represents  an  important  practical 
advance  of  NHA/OMT. 

The  preliminary  research  carried  out  at  Creative  Optics,  Inc.  during  Phase  I  indicates  the 
importance  of  numerical  condition  and  numerical  stability  both  for  NHA  and  to  a  lesser  extent 
for  NLLS  NHA.  The  choice  of  hyperbolic  functions  to  represent  radiances  is  well-motivated  by 
the  physics  of  the  atmosphere,  however,  it  is  also  well-known  that  not  all  energy  transport  in 
the  earth’s  atmosphere  is  radiative.  During  the  Phase  II  research  effort,  additional  basis 
functions  should  be  investigated  for  the  representation  of  radiance  data.  While  the  use  of 
hyperbolic  functions  as  a  part  of  the  NLLS  NHA  fitting  procedure  should  be  retained,  addition 
of  new  basis  functions  will  change  the  geometry  of  the  x2  surface.  This  will  yield 
improvements  in  the  numerical  condition  of  the  NLLS  NHA  and  improve  its  efficiency  for  the 
determination  of  global  minima  in  the  x2  surface  corresponding  to  physically  meaningful 
atmospheric  temperature  profiles. 
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Figure  4.  TOYS  Weight  FqiisJipn-  ^e  found  in  our  work  with  King’s  Differential  Inversion 

algorithm  that  atmospheric  weight  functions  are  well  represented  by  generalized  exponential 

weight  functions  shown  here  (solid  line)  for  m-0.742  and  P*  180.085  millitors,  plotted  together 

with  the  weight  function  for  TOYS  channel  16  (crosses). 
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5.  DATA  ANALYSIS  AND  DATA  EVALUATION 

5.1  Introduction.  During  the  Phase  1  research  effort  we  have  been  concerned  with  obtaining 
and/or  evaluating  suitable  radiance  data  sets  for  continued  development,  testing,  and  evaluation 
of  extensions  of  the  DI  algorithm  and  development  of  NHA/OMT.  This  establishes  a 
foundation  for  the  Phase  II  effort  which  will  be  directed  towards  developing  DI  and  NHA/OMT 
algorithms  to  a  state  suitable  for  the  routine  reduction  of  radiance  data  to  atmospheric 
temperatures  and  atmospheric  temperature  profiles. 

5.2  TOVS.  To  date,  data  from  the  Tiros  Operational  Vertical  Sounder  (TOVS)  has  been  the 
principal  resource  for  the  development  and  testing  of  DI  and  NHA/OMT  algorithms.  TOVs 
data  is  currently  available  in-house  at  Creative  Optics,  Inc.  with  15  scans  with  corresponding 
radiosonde  data  in  each  of  three  latitude  zones,  Zone  1  (polar,  latitude  =  60°  N  to  90°  N),  Zone 
2  (temperate,  latitude  *  30°  N  to  60°  N),  and  Zone  3  (equatorial,  latitude  =  0°  N  to  30°  N). 

A  summary  of  TOVS  instrumental  characteristics  is  given  in  Table  1.  Additional  TOVS 
data  is  readily  obtainable  for  support  of  the  Phase  II  research  effort.  [For  further  relevant 
information  on  TOVS,  see  (NOAA  1983;  Smith  et  al.  1979;  Weinreb  et  al.  1981)  and  references 
cited  therein.] 


TABLE  I.  Summary  of  TOVS  Instrumental  Characteristics. 


Half-power 


S  Channel  # 

Central  Wavenumberfcm'1) 

Plmb)  (zone  1) 

bandwid 

1 

668 

25.9 

3 

2 

679 

69.5 

10 

3 

691 

85. 

12 

4 

704 

400. 

16 

5 

716 

543. 

16 

6 

732 

780. 

16 

7 

748 

1000. 

16 

13 

2190 

992. 

23 

14 

2213 

663. 

23 

15 

2240 

400. 

23 

16 

2276 

172. 

23 

17 

2361 

25.9 

23 

-i\ 


5.3  HIS.  Data  were  obtained  from  the  High-resolution  Interferometer  Sounder  (HIS) 
developed  by  researchers  at  the  University  of  Wisconsin.  (See  Revercomb  et  al  1988.)  HIS  is  a 
Fourier-transform  spectrometer  operating  between  509  cm'1  and  2750  cm'1.  Relevant 
characteristics  of  the  HIS  instrument  are  given  in  Table  2.  A  set  of  15,000  radiances  and 
corresponding  weight  functions  have  been  obtained  from  HIS  and  a  significant  subset  of  them 
have  been  examined.  The  most  notable  qualitative  difference  between  the  HIS  data  and  the 
other  data  sets  we  have  examined  is  that  the  HIS  weight  functions  exhibit  multiple  maxima. 
(See  Figure  5a  for  a  representative  example.)  These  weight  functions  arise  due  to  the 
contributions  of  chemical  species  in  addition  to  carbon  dioxide.  DI  and  NHA/OMT  will  require 
some  generalization  in  order  to  reduce  radiance  data  such  as  these.  Until  such  theoretical 
generalization  is  accomplished,  reduction  of  HIS  data  will  be  restricted  to  those  wavelength 
intervals  where  only  absorption  of  a  single  species  is  important,  so  that  the  weight  function 
corresponding  to  that  wavelength  interval  exhibits  a  single  well-defined  maximum.  (See  Figure 
5b  for  an  example.)  Sufficient  wavelength  coverage  is  expected,  even  with  that  restriction,  to 
achieve  a  high  vertical  resolution  atmospheric  temperature  profile. 
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Figure  5.  HIS  Weighting  Functions.  Two  representative  weight  functions  for  the  High- 
resolution  Interferometer  Sounder  (HIS)  developed  by  the  University  of  Wisconsin.  The  weight 
function  in  (a)  exhibits  multiple  maxima  due  to  multiple  chemical  gaseous  species  contributing 
to  atmospheric  absorption  at  wavelength  frequency  of  731.209cm  .  The  single  peak  in  (b) 
corresponds  to  a  wavelength  interval  (centered  at  2244.391cm'1)  that  exhibits  a  single  well- 
defined  maximum. 
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TABLE  2.  Summar  of  HIS  Instrumental  Characteristics. 


HIS 

Band  I  590-1070  cm'1 
Band  II  1040-1930  cm'1 
Band  III  2070-2750  cm"1 


(Unapodized) 
Resolution 
.364  cm'1 
.641  cm'1 
.641  cm'1 


5.4  DMSP.  The  Defense  Meteorological  Satellite  Program  (DMSP)  satellites  carry  a  six- 
channel  microwave  radiometer  for  atmospheric  temperature  sounding.  The  relevant 
characteristics  of  this  instrument  are  listed  in  Table  3.  On  the  basis  of  the  instrumental 
characteristics  of  DMSP  and  published  weight  functions,  we  anticipate  that  reduction  of  DMSP 
by  DI  or  by  NHA/OMT  will  be  analogous  to  the  reduction  of  TOVS  data.  DMSP  data  could 
not  be  obtained  during  the  Phase  I  effort  because  radiances  are  only  available  from  Global 
Weather  Central  at  Offut  Air  Force  Base  in  a  36  bit  format  (as  required  for  management  of  the 
DMSP  database  by  Sperry-Univac  computers  used  at  Global  Weather  Central).  Data  can  be 
made  available  after  a  time  of  three  to  four  months  from  the  time  of  request,  which  is  the  time 
required  for  a  data  translation  program  to  be  made  by  workers  at  Global  Weather  Central. 


TABLE  3.  Summary  of  DMSP  Instrumental  Characteristics. 


DMSP 

Peaking 

Frequency 

Bandwidth 

NETD 

Calib 

Channel  # 

Heightfkm) 

(GHz)  . 

(MHz) 

(K) 

Uncertaintv(K) 

1 

0 

50.5 

400 

0.32 

2 

2 

53.2 

400 

0.24 

3 

6 

54.35 

400 

0.36 

0.24 

4 

10 

54.9 

400 

0.22 

0.13 

5 

30 

58.4 

115 

0.39 

0.06 

6 

16 

58.825 

400 

0.29 

0.24 

7 

22 

59.4 

250 

0.31 

0.13 
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A  number  of  conclusions  can  be  drawn  from  this  research  relevant  to  the  construction  of 
atmospheric  temperature  profiles  from  upwelling  radiance  data: 

6.1  Differential  Inversion  (Dll. 

•  When  the  series  representation  of  the  Planck  function  is  rapidly  convergent,  D1  can 
yield  multiple  atmospheric  temperature  values  suitable  for  the  construction  of  partial 
atmospheric  temperature  profiles. 

•  The  partial  temperature  profiles  determined  from  DI  are  suitable  for  bias-free 
determination  of  significant  meteorological  structures,  such  as  the  height  of  the  iropopause. 

•  DI  is  well-suited  to  radiance  data  sets  in  which  the  channel  weight  functions  sample 
the  atmosphere  finely  in  pressure. 

6.2  Optical  Measure  Theory  (OMT1  and  Nonlinear  Hyperbolic  Algorithm  fNHAL 

•  The  OMT  with  its  principal  algorithmic  element,  the  NHA,  have  been  implemented 
in  an  effective  FORTRAN  program  as  a  significant  step  towards  eventual  operational 
implementation  of  these  algorithms. 

•  The  algorithmic  properties  of  the  NHA  for  uniformly  noisy  radiance  data  and  for 
radiance  data  in  which  only  single  channels  are  corrupted  by  errors  have  been  characterized. 
The  ability  of  the  NHA  to  adaptively  process  radiance  data  with  single  corrupt  channels,  noted 
by  Leon  and  King,  has  been  confirmed  by  our  work. 

•  The  importance  for  NHA  of  uncertainties  in  pressure  values  at  which  atmospheric 
weight  functions  achieve  their  maximum  values  has  been  determined. 

•  Numerical  instability  of  the  NHA  is  now  well-characterized  and  understood  in  its 
relation  to  the  overall  stability  of  the  classic  problem  of  inversion  of  radiances  to  atmospheric 
temperatures. 


6.3  Nonlinear  Least  Squares  fNLLSl  Extension  of  NHA. 

•  A  nonlinear  least-squares  extension  of  the  NHA  has  been  demonstrated  and  shown  to 
be  an  appropriate  strategy  for  generating  as  many  hyperbolic  terms  in  the  NHA  as  the  data  will 
support. 

•  Improvement  of  the  NLLS  NHA  can  be  expected  as  a  result  of  application  of 
standard  algorithmic  techniques  appropriate  to  least-squares  curve-fitting  when  the  error  surface 
exhibits  multiple  secondary  minima. 

6.4  OMT  and  DI:  Bias-Free  Determination  of  Atmospheric  Temperature  Profiles. 

•  OMT  has  demonstrated  the  generation  of  a  bias-free  atmospheric  temperature  profile 
from  TOYS  data.  This  profile  exhibits  a  meteorological  character,  in  particular,  showing  a 
temperature  minimum  in  temperature  at  the  tropopause  altitude. 

•  The  algorithmic  basis  of  OMT,  from  transform  theory,  is  the  correct  approach  to  the 
calculation  of  atmospheric  temperature  profiles.  The  temperature  profiles  obtained  are  the 
smoothest  bias-free  temperatures  consistent  with  the  radiance  data. 


•  DI  and  OMT  are  complementary  bias-free  approaches  to  the  determination  of 
atmospheric  temperatures  from  upwelling  atmospheric  radiances.  (See  Figure  6.)  DI  calculates 
temperatures  at  some  set  of  pressure  values  in  the  atmosphere  with  no  underlying  assumption 
regarding  the  functional  form  of  the  radiance  data.  OMT  is  fundamentally  oriented  towards  the 
production  of  a  temperature  profile,  but  implicitly  assumes  that  the  dominant  mode  of  energy 
transport  in  the  atmosphere  is  radiative.  Additional  research  will  elucidate  the  relative  merits  of 
these  two  novel  approaches  to  remote  sensing  of  atmospheric  temperatures. 
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BIAS-FREE  INVERSE  TRANSFORMATION- 
DEMONSTRATION  ATMOSPHERIC  TEMPERATURE  PROFILE  USING 
OPTICAL  MEASURE  THEORY  APPLIED  TO  TOVS  DATA 
(NONLINEAR  LEAST  SQUARES  IMPLEMENTATION) 


Pressure(mb) 


220  225  230  235  240  245  250  255  260  265  270 
Temperature(K) 


Figure  6.  Comparison  of  Two  Methods  of  Bias-Free  Temperature  Inversion.  Shown  is  a 
comparison  of  temperature  profiling  methods  implemented  by  Creative  Optics,  Inc  — 
Differential  Inversion  (DI)  and  Optical  Measure  Theory  (OMT)--together  with  corresponding 
RAOB  data.  DI  and  OMT  represent  complementary  bias-free  methods  for  the  determination  of 
atmospheric  temperatures.  DI  yields  atmospheric  temperatures  with  no  assumption  as  to  an 
underlying  functional  form,  while  OMT  produces  the  smoothest  temperature  profile  consistent 
with  the  radiance  data.  In  conjunction,  these  two  methods  yield  important  new  information  on 
atmospheric  temperatures  from  upweliing  radiance  data. 
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APPENDIX  A 

THREE-POINT  TOVS  TEMPERATURE  INFERENCE  USING  DIFFERENTIAL  INVERSION 

We  were  able  to  determine  temperatures  from  the  TOVS  data  at  three  levels  in  the  atmosphere 
in  the  4.3jim  wavelength  band.  The  values  of  the  pressure  at  these  levels  are  172.3,  400,  663 
millibars.  The  figures  are  plots  of  the  values  for  temperature  from  differential  inversion 
together  with  RAOB  measurements  against  pressure.  The  RAOB  measurements  are  given  by  the 
solid  line.  Note  that  point  on  both  sides  of  the  tropopause  were  obtained 
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